#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "IDoubleMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class DoubleSymmetricMatrix : public IDoubleMatrix {
	internal:
		double * a;
	private:
		long * ipiv;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;

		DoubleSymmetricMatrix(DoubleSymmetricMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;

			int size = m_rownum * (m_rownum + 1) / 2;
			a = new double[size];

			for (int i = 0; i < size; i++) {
				a[i] = orig->a[i];
			}
			ipiv = new long[m_rownum];

			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}

			m_isFactorized = orig->m_isFactorized;

		};
	public:
		DoubleSymmetricMatrix(int rownum, int colnum) {
			if (rownum != colnum) {
				throw gcnew LatoolException("Number of rows and columns must be equal for symmetric matrix.");
			}
			m_rownum = rownum;
			m_colnum = colnum;

			int size = rownum * (rownum + 1) / 2;
			a = new double[size];

			for (int i = 0; i < size; i++) {
				a[i] = 0.0;
			}
      ipiv = new long[rownum];

			m_isFactorized = false;
		};

		~DoubleSymmetricMatrix(){
			delete[] a;
		};

		virtual IDoubleMatrix ^ Clone() {		
			return (IDoubleMatrix^) gcnew DoubleSymmetricMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property double default[int, int] {
			double get(int i, int j) {
				if (i > j) {
					int temp = i;
					i = j;
					j = temp;
				}
				return a[i + j * (j + 1) / 2];
			}
			void set(int i, int j, double value) {
				if (i > j) {
					int temp = i;
					i = j;
					j = temp;
				}
				a[i + j * (j + 1) / 2] = value;
				m_isFactorized = false;
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::DoubleSymmetric; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		}

		virtual void Invert() {

			char uplo = 'U';
			long n = m_rownum;
			long info;
				
			double * work = new double[Math::Max(1, n)];

			dsptri_(&uplo, &n, a, ipiv, work, &info);

			delete[] work;

			m_isFactorized = false;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("the i-th diagonal element of the Cholesky factor (and therefore the factor itself) is zero, and the inversion could not be completed.");
			}

		};

		virtual void Factorize() {

			char uplo = 'U';
			long n = m_rownum;
			long info;

			dsptrf_(&uplo, &n, a, ipiv, &info);

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("The factorization has been completed, but D is exactly singular. Division by 0 will occur if you use D for solving a system of linear equations.");
			}

			m_isFactorized = true;

		};

		virtual void Solve(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			char uplo = 'U';
			long n = m_rownum;
			long nrhs = b->ColNum;
			long ldb = Math::Max(1, n);
			long info;

			dsptrs_(&uplo, &n, &nrhs, a, ipiv, gb->dat, &ldb, &info);

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			char uplo = 'U';
			long n = m_rownum;
			long nrhs = b->ColNum;
			long ldb = Math::Max(1, n);
			long info;

			dspsv_(&uplo, &n, &nrhs, a, ipiv, gb->dat, &ldb, &info);

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("The factorization has been completed, but D is exactly singular, so the solution could not be computed.");
			}

			m_isFactorized = true;

		};

		void SolveSymmetricEigenproblem(array<double>^ values, IMatrix^ vectors) {

			char jobz = 'V';
			char uplo = 'U';
			long n = m_rownum;
			double * work = new double[3 * n];
			long ldz = n;
			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;
			double * z = ((DoubleGeneralMatrix ^) vectors)->dat;
			long info = 0;
		
			dspev_(&jobz, &uplo, &n, a, w, z, &ldz, work, &info);

			delete[] work;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException(String::Format("{i} elements of an intermediate tridiagonal form did not converge to zero.", info));
			}

		};

		void SolveSymmetricEigenproblem(array<double>^ values) {

			char jobz = 'N';
			char uplo = 'U';
			long n = m_rownum;
			double * work = new double[3 * n];
			long ldz = n;
			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;
			double z = 0;
			long info = 0;
		
			dspev_(&jobz, &uplo, &n, a, w, &z, &ldz, work, &info);

			delete[] work;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException(String::Format("{i} elements of an intermediate tridiagonal form did not converge to zero.", info));
			}

		};

		void SolveGeneralizedSymmetricDefiniteEigenproblem(IMatrix^ b, array<double>^ values, IMatrix^ vectors) {
			
			long itype = 1;
			char jobz = 'V';
			char uplo = 'U';
			long n = m_rownum;
			double * ap = a;
			double * bp = ((DoubleSymmetricMatrix ^) b)->a;

			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;

			double * z = ((DoubleGeneralMatrix ^) vectors)->dat;
			long ldz = n;
			double * work = new double[3 * n];
			long info = 0;

			dspgv_(&itype, &jobz, &uplo, &n, ap, bp, w, z, &ldz, work, &info);

			delete[] work;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("off-diagonal elements of an intermediate tridiagonal did not converge to zero");
			}
		};

		void SolveGeneralizedSymmetricDefiniteEigenproblem(IMatrix^ b, array<double>^ values) {
			
			long itype = 1;
			char jobz = 'N';
			char uplo = 'U';
			long n = m_rownum;
			double * ap = a;
			double * bp = ((DoubleSymmetricMatrix ^) b)->a;

			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;

			double z;
			long ldz = n;
			double * work = new double[3 * n];
			long info = 0;

			dspgv_(&itype, &jobz, &uplo, &n, ap, bp, w, &z, &ldz, work, &info);

			delete[] work;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("off-diagonal elements of an intermediate tridiagonal did not converge to zero");
			}

		};

		void SPEV(IMatrix ^ eigenValue, IMatrix ^ eigenVector, long * info) {
			char jobz = 'V';
			char uplo = 'U';
			long n = m_rownum;
			double * work = new double[3 * n];
			long ldz = n;
			double * w = ((DoubleGeneralMatrix ^) eigenValue)->dat;
			double * z = ((DoubleGeneralMatrix ^) eigenVector)->dat;
		
			dspev_(&jobz, &uplo, &n, a, w, z, &ldz, work, info);

			delete[] work;


		};

		void SPGV(IMatrix ^ b, IMatrix ^ eigenValue, IMatrix ^ eigenVector, long * info) {
			long itype = 1;
			char jobz = 'V';
			char uplo = 'U';
			long n = m_rownum;
			double * ap = a;
			double * bp = ((DoubleSymmetricMatrix ^) b)->a;
			double * w = ((DoubleGeneralMatrix ^) eigenValue)->dat;
			double * z = ((DoubleGeneralMatrix ^) eigenVector)->dat;
			long ldz = n;
			double * work = new double[3 * n];

			dspgv_(&itype, &jobz, &uplo, &n, ap, bp, w, z, &ldz, work, info);

			delete[] work;
		};


	};
}